#include "boundConds.h"
#include <iostream>

using namespace std;

BoundConds::BoundConds(const Config &cfg, const Geometry &geo, Fluid &flow)
    : equsType{cfg.GetEquationType()},
      Bdescrp{geo.BoundDescrp}, si{geo.si}, sj{geo.sj},
      cv{flow.GetCV()}, dvE{flow.GetdvEuler()},
      dvNS{flow.GetdvNS()}, press{flow.GetPressure()},
      ua{flow.GetuInf()}, va{flow.GetvInf()}, pa{flow.GetpInf()},
      rhoa{flow.GetRhoInf()}, mach{cfg.GetMachNumber()},
      il{geo.il}, jl{geo.jl}, i2{geo.i2}, j2{geo.j2},
      imax{il + 1}, jmax{jl + 1},
      cutExchgCV(i2, j2), cutExchgDVE(i2, j2),
      calcDvOneEu(flow.GetRgas()),
      calcDvOneNS(flow.GetRgas(), flow.GetSutherlandMuref(),
                  flow.GetSutherlandTref(), flow.GetSutherlandTs())

{
}

void BoundConds::BoundaryConditionsEuler(void)
{
  // fill dummy cells at the corners: (1,1), (il,1), (1,jl), (il,jl)
  ExtraAtCorner(cv(1, 1), cv(1, 2), cv(2, 1));
  ExtraAtCorner(cv(il, 1), cv(il, 2), cv(i2, 1));
  ExtraAtCorner(cv(1, jl), cv(1, j2), cv(2, jl));
  ExtraAtCorner(cv(il, jl), cv(il, j2), cv(i2, jl));
  calcDvOneEu(cv(1, 1), dvE(1, 1), press(1, 1));
  calcDvOneEu(cv(il, 1), dvE(il, 1), press(il, 1));
  calcDvOneEu(cv(1, jl), dvE(1, jl), press(1, jl));
  calcDvOneEu(cv(il, jl), dvE(il, jl), press(il, jl));

  // loop over all boundaries with the exception of walls and symmetry planes
  int nbound = Bdescrp.size();
  for (int i = 0; i < nbound; ++i)
  {
    switch (Bdescrp[i].BType)
    {
    case Bound_t::Farfield:
      // printArray(cout, cv, "before farfield bound   cv");
      // 没有外拓  依赖变量取的都是内部的
      Farfield(Bdescrp[i], dvE);
      // printArray(cout, cv, "after farfield bound set   cv");
      break;
    case Bound_t::InnerBound:
      cutExchgCV(Bdescrp[i], cv);
      cutExchgDVE(Bdescrp[i], dvE);
      break;
    default:
      // cerr << "this part only support Farfield cut now\n"
      //      << endl;
      // exit(EXIT_FAILURE);
      break;
    }
  }

  // solid walls and symmetry (treated last because they should dominate)
  for (int i = 0; i < nbound; ++i)
  {
    switch (Bdescrp[i].BType)
    {
    case Bound_t::Wall:
      // 外拓一层  没有使用依赖变量
      WallEuler(Bdescrp[i]);
      break;
    case Bound_t::Symmetry:
      // 外拓一层 没有使用依赖变量
      Symmetry(Bdescrp[i]);
      break;
    default:
      // cerr << "this part only support Farfield cut now\n"
      //      << endl;
      // exit(EXIT_FAILURE);
      break;
    }
  }
  DependentVarsBoundaryEu();
}

void BoundConds::BoundaryConditionsNS(void)
{
  // fill dummy cells at the corners: (1,1), (il,1), (1,jl), (il,jl)
  ExtraAtCorner(cv(1, 1), cv(1, 2), cv(2, 1));
  ExtraAtCorner(cv(il, 1), cv(il, 2), cv(i2, 1));
  ExtraAtCorner(cv(1, jl), cv(1, j2), cv(2, jl));
  ExtraAtCorner(cv(il, jl), cv(il, j2), cv(i2, jl));
  calcDvOneNS(cv(1, 1), dvNS(1, 1), press(1, 1));
  calcDvOneNS(cv(il, 1), dvNS(il, 1), press(il, 1));
  calcDvOneNS(cv(1, jl), dvNS(1, jl), press(1, jl));
  calcDvOneNS(cv(il, jl), dvNS(il, jl), press(il, jl));

  // loop over all boundaries with the exception of walls and symmetry planes
  int nbound = Bdescrp.size();
  for (int i = 0; i < nbound; ++i)
  {
    switch (Bdescrp[i].BType)
    {
    case Bound_t::Farfield:
      // printArray(cout, cv, "before farfield bound   cv");
      // 没有外拓  依赖变量取的都是内部的
      Farfield(Bdescrp[i], dvNS);
      // printArray(cout, cv, "after farfield bound set   cv");
      break;
    case Bound_t::InnerBound:
      cutExchgCV(Bdescrp[i], cv);
      cutExchgDVNS(Bdescrp[i], dvNS);
      break;
    default:
      // cerr << "this part only support Farfield cut now\n"
      //      << endl;
      // exit(EXIT_FAILURE);
      break;
    }
  }

  // solid walls and symmetry (treated last because they should dominate)
  for (int i = 0; i < nbound; ++i)
  {
    switch (Bdescrp[i].BType)
    {
    case Bound_t::Wall:
      // 外拓一层  没有使用依赖变量
      WallNS(Bdescrp[i]);
      break;
    case Bound_t::Symmetry:
      // 外拓一层 没有使用依赖变量
      Symmetry(Bdescrp[i]);
      break;
    default:
      // cerr << "this part only support Farfield symmetry now\n"
      //      << endl;
      // exit(EXIT_FAILURE);
      break;
    }
  }
  DependentVarsBoundaryNS();
}

void BoundConds::WallNS(const BoundDescrp_t &Bdescrp)
{
  // const int(&trans)[4] = Bdescrp.isIface ? transI : transJ;
  Transij trans(Bdescrp.isIface);
  const int dum1{Bdescrp.dum1},
      dum2{Bdescrp.dum2},
      ins1{Bdescrp.ins1},
      ins2{Bdescrp.ins2};
  int beg{Bdescrp.beg},
      end{Bdescrp.end};

  // extend to dummy cells (1st layer)
  beg = (beg == 2) ? 1 : beg;
  end = (end == i2 || end == j2) ? end + 1 : end;

  int idum1, jdum1, idum2, jdum2, iins1, jins1, iins2, jins2;
  for (int i = beg; i <= end; ++i)
  {
    trans(i, dum1, idum1, jdum1);
    trans(i, ins1, iins1, jins1);
    trans(i, dum2, idum2, jdum2);
    trans(i, ins2, iins2, jins2);

    ConsVar_t &cvDum1{cv(idum1, jdum1)}, &cvDum2{cv(idum2, jdum2)};
    ConsVar_t &cvIns1{cv(iins1, jins1)}, &cvIns2{cv(iins2, jins2)};

    cvDum1.rho = cvIns1.rho;
    cvDum1.rhou = -cvIns1.rhou;
    cvDum1.rhov = -cvIns1.rhov;
    cvDum1.rhoE = cvIns1.rhoE;

    cvDum2.rho = cvIns2.rho;
    cvDum2.rhou = -cvIns2.rhou;
    cvDum2.rhov = -cvIns2.rhov;
    cvDum2.rhoE = cvIns2.rhoE;

    // DependentVarsOne(cvDum1, dv(idum1, jdum1));
    // DependentVarsOne(cvDum2, dv(idum2, jdum2));
  }
}

void BoundConds::WallEuler(const BoundDescrp_t &Bdescrp)
{
  Transij trans(Bdescrp.isIface);
  const int dum1{Bdescrp.dum1}, dum2{Bdescrp.dum2},
      ins1{Bdescrp.ins1}, ins2{Bdescrp.ins2};
  int beg{Bdescrp.beg}, end{Bdescrp.end};

  // extend to dummy cells (1st layer)
  beg = (beg == 2) ? 1 : beg;
  end = (end == i2 || end == j2) ? end + 1 : end;

  int idum1, jdum1, idum2, jdum2, iins1, jins1, iins2, jins2;
  for (int i = beg; i <= end; ++i)
  {
    trans(i, dum1, idum1, jdum1);
    trans(i, ins1, iins1, jins1);
    trans(i, dum2, idum2, jdum2);
    trans(i, ins2, iins2, jins2);

    ConsVar_t &cvDum1{cv(idum1, jdum1)}, &cvDum2{cv(idum2, jdum2)};
    ConsVar_t &cvIns1{cv(iins1, jins1)}, &cvIns2{cv(iins2, jins2)};

    cvDum1.rho = 2.0 * cvIns1.rho - cvIns2.rho;
    cvDum1.rhou = 2.0 * cvIns1.rhou - cvIns2.rhou;
    cvDum1.rhov = 2.0 * cvIns1.rhov - cvIns2.rhov;
    cvDum1.rhoE = 2.0 * cvIns1.rhoE - cvIns2.rhoE;

    cvDum2.rho = 2.0 * cvDum1.rho - cvIns1.rho;
    cvDum2.rhou = 2.0 * cvDum1.rhou - cvIns1.rhou;
    cvDum2.rhov = 2.0 * cvDum1.rhov - cvIns1.rhov;
    cvDum2.rhoE = 2.0 * cvDum1.rhoE - cvIns1.rhoE;

    // DependentVarsOne(cvDum1, dv(idum1, jdum1));
    // DependentVarsOne(cvDum2, dv(idum2, jdum2));
  }
}

void BoundConds::Symmetry(const BoundDescrp_t &Bdescrp)
{
  const bool isIface{Bdescrp.isIface};
  Transij trans(isIface);
  const int dum1{Bdescrp.dum1}, dum2{Bdescrp.dum2}, ins1{Bdescrp.ins1},
      ins2{Bdescrp.ins2};
  int beg{Bdescrp.beg}, end{Bdescrp.end};

  const VarArray<Area_t, 2, 0> &S{isIface ? si : sj};
  int indBface{(Bdescrp.Bn < 2) ? ins1 : dum1};

  int iBface, jBface;
  int idum1, jdum1, idum2, jdum2, iins1, jins1, iins2, jins2;
  real_t qn, sxn, syn, ds;
  for (int i = beg; i <= end; ++i)
  {
    trans(i, indBface, iBface, jBface);
    trans(i, dum1, idum1, jdum1);
    trans(i, ins1, iins1, jins1);
    trans(i, dum2, idum2, jdum2);
    trans(i, ins2, iins2, jins2);

    const Area_t &bs{S(iBface, jBface)};
    // 外拓的话就会出问题   面积在最边界有可能为0
    ds = sqrt(bs.sx * bs.sx + bs.sy * bs.sy);
    sxn = bs.sx / ds;
    syn = bs.sy / ds;

    ConsVar_t &cvDum1{cv(idum1, jdum1)}, &cvDum2{cv(idum2, jdum2)};
    const ConsVar_t &cvIns1{cv(iins1, jins1)}, &cvIns2{cv(iins2, jins2)};
    // 垂直面的速度
    qn = cvIns1.rhou * sxn + cvIns1.rhov * syn;

    cvDum1.rho = cvIns1.rho;
    cvDum1.rhou = cvIns1.rhou - 2.0 * qn * sxn;
    cvDum1.rhov = cvIns1.rhov - 2.0 * qn * syn;
    cvDum1.rhoE = cvIns1.rhoE;

    qn = cvIns2.rhou * sxn + cvIns2.rhov * syn;

    cvDum2.rho = cvIns2.rho;
    cvDum2.rhou = cvIns2.rhou - 2.0 * qn * sxn;
    cvDum2.rhov = cvIns2.rhov - 2.0 * qn * syn;
    cvDum2.rhoE = cvIns2.rhoE;

    // DependentVarsOne(cvDum1, dv(idum1, jdum1));
    // DependentVarsOne(cvDum2, dv(idum2, jdum2));
  }

  // 下面为书上建议外拓的 一层
  // /**
  if (beg == 2)
  {
    // 取内边界的面  直接平移过去
    trans(2, indBface, iBface, jBface);
    trans(1, dum1, idum1, jdum1);
    trans(1, ins1, iins1, jins1);
    trans(1, dum2, idum2, jdum2);
    trans(1, ins2, iins2, jins2);

    const Area_t &bs{S(iBface, jBface)};
    // 外拓的话就会出问题
    ds = sqrt(bs.sx * bs.sx + bs.sy * bs.sy);
    sxn = bs.sx / ds;
    syn = bs.sy / ds;

    ConsVar_t &cvDum1{cv(idum1, jdum1)}, &cvDum2{cv(idum2, jdum2)};
    const ConsVar_t &cvIns1{cv(iins1, jins1)}, &cvIns2{cv(iins2, jins2)};

    qn = cvIns1.rhou * sxn + cvIns1.rhov * syn;

    cvDum1.rho = cvIns1.rho;
    cvDum1.rhou = cvIns1.rhou - 2.0 * qn * sxn;
    cvDum1.rhov = cvIns1.rhov - 2.0 * qn * syn;
    cvDum1.rhoE = cvIns1.rhoE;

    qn = cvIns2.rhou * sxn + cvIns2.rhov * syn;

    cvDum2.rho = cvIns2.rho;
    cvDum2.rhou = cvIns2.rhou - 2.0 * qn * sxn;
    cvDum2.rhov = cvIns2.rhov - 2.0 * qn * syn;
    cvDum2.rhoE = cvIns2.rhoE;

    // DependentVarsOne(cvDum1, dv(idum1, jdum1));
    // DependentVarsOne(cvDum2, dv(idum2, jdum2));
  }
  // **/
  if ((end == i2) || (end == j2))
  {
    // 取内边界的面  直接平移过去
    trans(end, indBface, iBface, jBface);
    trans(end + 1, dum1, idum1, jdum1);
    trans(end + 1, ins1, iins1, jins1);
    trans(end + 1, dum2, idum2, jdum2);
    trans(end + 1, ins2, iins2, jins2);

    const Area_t &bs{S(iBface, jBface)};
    // 外拓的话就会出问题
    ds = sqrt(bs.sx * bs.sx + bs.sy * bs.sy);
    sxn = bs.sx / ds;
    syn = bs.sy / ds;

    ConsVar_t &cvDum1{cv(idum1, jdum1)}, &cvDum2{cv(idum2, jdum2)};
    const ConsVar_t &cvIns1{cv(iins1, jins1)}, &cvIns2{cv(iins2, jins2)};

    qn = cvIns1.rhou * sxn + cvIns1.rhov * syn;

    cvDum1.rho = cvIns1.rho;
    cvDum1.rhou = cvIns1.rhou - 2.0 * qn * sxn;
    cvDum1.rhov = cvIns1.rhov - 2.0 * qn * syn;
    cvDum1.rhoE = cvIns1.rhoE;

    qn = cvIns2.rhou * sxn + cvIns2.rhov * syn;

    cvDum2.rho = cvIns2.rho;
    cvDum2.rhou = cvIns2.rhou - 2.0 * qn * sxn;
    cvDum2.rhov = cvIns2.rhov - 2.0 * qn * syn;
    cvDum2.rhoE = cvIns2.rhoE;

    // DependentVarsOne(cvDum1, dv(idum1, jdum1));
    // DependentVarsOne(cvDum2, dv(idum2, jdum2));
  }
}

template <typename T>
void BoundConds::Farfield(const BoundDescrp_t &Bdescrp, const VarArray<T, 2, 0> &dv)
{
  const int dum1{Bdescrp.dum1}, dum2{Bdescrp.dum2},
      ins1{Bdescrp.ins1},
      beg{Bdescrp.beg}, end{Bdescrp.end},
      indBface{(Bdescrp.Bn < 2) ? ins1 : dum1};

  Transij trans(Bdescrp.isIface);
  const VarArray<Area_t, 2, 0> &S{Bdescrp.isIface ? si : sj};
  real_t ds, sxn, syn, sgn{(Bdescrp.Bn < 2) ? 1.0 : -1.0};

  int i1, j1, iins, jins, idum1, jdum1, idum2, jdum2;
  const ConsVar_t cvInf{rhoa, rhoa * ua, rhoa * va,
                        pa / gam1 + 0.5 * rhoa * (ua * ua + va * va)};
  real_t pd, ud, vd, rhod, cd, qn, crho0;
  real_t pb;

  for (int i = beg; i <= end; ++i)
  {
    trans(i, indBface, i1, j1);
    // 注意面与单元索引的不同
    // 使面指向外部  为了区分内流和外流 inflow/outflow
    const Area_t &bs{S(i1, j1)};
    ds = sqrt(bs.sx * bs.sx + bs.sy * bs.sy);
    sxn = sgn * bs.sx / ds;
    syn = sgn * bs.sy / ds;

    trans(i, ins1, iins, jins);
    const ConsVar_t &cvIns{cv(iins, jins)};
    // 流场内部的信息
    rhod = cvIns.rho;
    ud = cvIns.rhou / rhod;
    vd = cvIns.rhov / rhod;
    // qd = ud*ud + vd*vd;
    pd = press(iins, jins);
    cd = dv(iins, jins).c;
    crho0 = cd * rhod;
    qn = sxn * ud + syn * vd;

    trans(i, dum1, idum1, jdum1);
    ConsVar_t &cvDum1{cv(idum1, jdum1)};
    // 亚音速流
    if (mach < 1.0)
    {
      // inflow
      if (qn < 0.0)
      {
        pb = 0.5 * (pa + pd - crho0 * (sxn * (ua - ud) + syn * (va - vd)));
        cvDum1.rho = rhoa + (pb - pa) / (cd * cd);
        cvDum1.rhou = cvDum1.rho * (ua - sxn * (pa - pb) / crho0);
        cvDum1.rhov = cvDum1.rho * (va - syn * (pa - pb) / crho0);
        cvDum1.rhoE = pb / gam1 + 0.5 * (cvDum1.rhou * cvDum1.rhou + cvDum1.rhov * cvDum1.rhov) / cvDum1.rho;
      }
      // outflow
      else
      {
        pb = pa;
        cvDum1.rho = rhod + (pb - pd) / (cd * cd);
        cvDum1.rhou = cvDum1.rho * (ud + sxn * (pd - pb) / crho0);
        cvDum1.rhov = cvDum1.rho * (vd + syn * (pd - pb) / crho0);
        cvDum1.rhoE = pb / gam1 + 0.5 * (cvDum1.rhou * cvDum1.rhou + cvDum1.rhov * cvDum1.rhov) / cvDum1.rho;
      }
    }
    // 超音速流
    else
    {
      cvDum1 = (qn < 0.0) ? cvInf : cvIns;
    }
    // 第二层虚单元
    trans(i, dum2, idum2, jdum2);
    ConsVar_t &cvDum2{cv(idum2, jdum2)};
    cvDum2.rho = 2.0 * cvDum1.rho - cvIns.rho;
    cvDum2.rhou = 2.0 * cvDum1.rhou - cvIns.rhou;
    cvDum2.rhov = 2.0 * cvDum1.rhov - cvIns.rhov;
    cvDum2.rhoE = 2.0 * cvDum1.rhoE - cvIns.rhoE;

    // DependentVarsOne(cvDum1, dv(idum1, jdum1));
    // DependentVarsOne(cvDum2, dv(idum2, jdum2));
  }
}

void BoundConds::DependentVarsBoundaryEu(void)
{
  for (int j = 0; j <= jmax; ++j)
  {
    calcDvOneEu(cv(0, j), dvE(0, j), press(0, j));
    calcDvOneEu(cv(1, j), dvE(1, j), press(1, j));
    calcDvOneEu(cv(il, j), dvE(il, j), press(il, j));
    calcDvOneEu(cv(imax, j), dvE(imax, j), press(imax, j));
  }
  for (int i = 2; i <= i2; ++i)
  {
    calcDvOneEu(cv(i, 0), dvE(i, 0), press(i, 0));
    calcDvOneEu(cv(i, 1), dvE(i, 1), press(i, 1));
    calcDvOneEu(cv(i, jl), dvE(i, jl), press(i, jl));
    calcDvOneEu(cv(i, jmax), dvE(i, jmax), press(i, jmax));
  }
}

void BoundConds::DependentVarsBoundaryNS(void)
{
  for (int j = 0; j <= jmax; ++j)
  {
    calcDvOneNS(cv(0, j), dvNS(0, j), press(0, j));
    calcDvOneNS(cv(1, j), dvNS(1, j), press(1, j));
    calcDvOneNS(cv(il, j), dvNS(il, j), press(il, j));
    calcDvOneNS(cv(imax, j), dvNS(imax, j), press(imax, j));
  }
  for (int i = 2; i <= i2; ++i)
  {
    calcDvOneNS(cv(i, 0), dvNS(i, 0), press(i, 0));
    calcDvOneNS(cv(i, 1), dvNS(i, 1), press(i, 1));
    calcDvOneNS(cv(i, jl), dvNS(i, jl), press(i, jl));
    calcDvOneNS(cv(i, jmax), dvNS(i, jmax), press(i, jmax));
  }
}